/*
 * Test_Compare_EquilibriumSolver_TotalEquilibriumSolver.h
 *
 *  Created on: 10 Nov 2011
 *      Author: allan
 */

#ifndef TEST_COMPARE_EQUILIBRIUMSOLVER_TOTALEQUILIBRIUMSOLVER_H_
#define TEST_COMPARE_EQUILIBRIUMSOLVER_TOTALEQUILIBRIUMSOLVER_H_

#include <GeoReact.h>
#include <Eigen/Core>

#include <boost/timer.hpp>

void Test_EquilibriumSolver();
void Test_TotalEquilibriumSolver();
	
void Test_Compare_EquilibriumSolver_TotalEquilibriumSolver()
{
	Test_EquilibriumSolver();
	Test_TotalEquilibriumSolver();
}

void Test_EquilibriumSolver()
{
	GeoManager manager;
	
	manager.SetAqueousSpecies({"CO2(aq)", "CO3[2-]", "CaCO3(aq)", "CaCl2(aq)", "CaCl[+]", 
		                       "CaHCO3[+]", "CaOH[+]", "Ca[2+]", "Cl[-]", "H2O", "HCO3[-]", 
		                       "HCl(aq)", "H[+]", "NaCO3[-]", "NaCl(aq)", "NaHCO3(aq)", "NaOH(aq)", "Na[+]", "OH[-]"});
	
	manager.SetGaseousSpecies({"H2O(g)", "CO2(g)"});
	
	GeoSystem system = manager.AssembleGeoSystem();
	
//	cout << system << endl;
	
	GeoEquilibrator equilibrator(system);
	
	equilibrator.ImposeElementMass("Ca[2+]",  4.9653284092e-04);
	equilibrator.ImposeElementMass("Cl[-]",   6.0000004081e+00);
	equilibrator.ImposeElementMass("H2O",     5.3996455542e+01);
	equilibrator.ImposeElementMass("HCO3[-]", 1.5295649624e+00);
	equilibrator.ImposeElementMass("H[+]",    1.5285718966e+00);
	equilibrator.ImposeElementMass("Na[+]",   6.0000004082e+00);
	
	GeoState state(system);
	
	state.SetTemperature(60.0);
	state.SetPressure(80.0);
	
	VectorXd n(21);
	
	n << 3.9153286345e-01,	     5.0296390890e-10,	     5.8557224907e-12,	     6.9518264236e-05,	     4.1466759670e-05,	     3.3866823685e-07,	     6.5492485554e-15,
	     1.9884210456e-04,	     3.9724158098e+00,	     5.5507965823e+01,	     7.0002104892e-04,	     2.5381645034e-04,	     8.5484607369e-04,	     3.5776285531e-10,
	     2.0271502786e+00,	     1.0286326672e-03,	     4.1426411868e-11,	     3.9718214965e+00,	     7.0689859803e-11,	     1.6139321634e-02,	     1.1361167387e+00;

	const VectorXd nCopy = n;
	
	equilibrator.Equilibrate(state, n);
	
//	cout << state << endl;
	
	cout << endl << "Iterations(EquilibriumSolver): " << equilibrator.GetDiagnostic().iterations << endl;
	
	boost::timer timer;
	
	for(unsigned i = 0; i < 100; ++i)
	{
		n = nCopy;
		equilibrator.Equilibrate(state, n);
	}
	
	cout << "ElapsedTime(EquilibriumSolver): " << timer.elapsed() << endl;
}

void Test_TotalEquilibriumSolver()
{
	double T = 60.0;
	double P = 80.0;
	
	GeoManager manager;
	
	manager.SetAqueousSpecies({"CO2(aq)", "CO3[2-]", "CaCO3(aq)", "CaCl2(aq)", "CaCl[+]", 
		                       "CaHCO3[+]", "CaOH[+]", "Ca[2+]", "Cl[-]", "H2O", "HCO3[-]", 
		                       "HCl(aq)", "H[+]", "NaCO3[-]", "NaCl(aq)", "NaHCO3(aq)", "NaOH(aq)", "Na[+]", "OH[-]"});
	
	manager.SetGaseousSpecies({"H2O(g)", "CO2(g)"});
	
	GeoSystem system = manager.AssembleGeoSystem();
	
	vector<string> jSpecies = {"Ca[2+]", "Cl[-]", "H2O", "HCO3[-]", "H[+]", "Na[+]"};
	vector<string> eSpecies = {"CO2(aq)", "CO3[2-]", "CaCO3(aq)", "CaCl2(aq)", "CaCl[+]", "CaHCO3[+]", "CaOH[+]", "HCl(aq)", "NaCO3[-]", "NaCl(aq)", "NaHCO3(aq)", "NaOH(aq)", "OH[-]", "H2O(g)", "CO2(g)"};

	TotalEquilibriumSolver solver(system.GetMultiphase(), system.GetEquilibriumReactions(), jSpecies, eSpecies);
	
	VectorXd uj(6);
	
	uj << 4.9653284092e-04,          6.0000004081e+00,          5.3996455542e+01,          1.5295649624e+00,          1.5285718966e+00,          6.0000004082e+00;
	
	VectorXd n(21);
	
	n << 3.9153286345e-01,	     5.0296390890e-10,	     5.8557224907e-12,	     6.9518264236e-05,	     4.1466759670e-05,	     3.3866823685e-07,	     6.5492485554e-15,
	     1.9884210456e-04,	     3.9724158098e+00,	     5.5507965823e+01,	     7.0002104892e-04,	     2.5381645034e-04,	     8.5484607369e-04,	     3.5776285531e-10,
	     2.0271502786e+00,	     1.0286326672e-03,	     4.1426411868e-11,	     3.9718214965e+00,	     7.0689859803e-11,	     1.6139321634e-02,	     1.1361167387e+00;
	
	const VectorXd nCopy = n;
	
	auto diagnostic = solver.Solve(uj, T, P, n);

	cout << "Iterations(TotalEquilibriumSolver): " << diagnostic.iterations << endl;
	
	boost::timer timer;
	
	for(unsigned i = 0; i < 100; ++i)
	{
		n = nCopy;
		solver.Solve(uj, T, P, n);
	}
	
	cout << "ElapsedTime(TotalEquilibriumSolver): " << timer.elapsed() << endl;
}




#endif /* TEST_COMPARE_EQUILIBRIUMSOLVER_TOTALEQUILIBRIUMSOLVER_H_ */
